Single-cell landscape of immunocytes in patients with extrahepatic cholangiocarcinoma

Background The intricate landscape of immunocytes in the tumor microenvironment (TME) is fundamental to immunotherapy but notably under-researched in extrahepatic cholangiocarcinoma (ECCA). Methods Single-cell RNA sequencing technology was conducted to make an in-depth analysis of immunocytes from matched tumor tissues, paratumor tissues and peripheral blood from ECCA patients. The potential cellular interactions between two cell populations were analyzed with software CellPhoneDB (v2.1.7). Results We obtained 13526 cells and characterized the transcriptomes and heterogeneity of different clusters and subclusters of immunocytes from ECCA, including CD4+ T cells, CD8+ T cells, B cells and myeloid immunocytes. We observed the rarely described immunocyte subclusters "intermediate" exhausted CD8+ T (CD8+ Tex) cells and “nonclassic” plasmacytes (CD27+ CD138+ CD38−). In addition, we identified potential immunotherapy targets, for example, ACP5, MAGEH1, TNFRSF9 and CCR8 for Tregs and MT1 for CD8+ Tex cells. We also found strong cellular interactions among Treg cells, M2 macrophages and CD8+ Tex cells through ligand–receptor analysis, implying that potential cellular cross-linkage promoted the immunosuppressive nature of the TME. Conclusions In a word, our study illuminated the components of the TME and revealed potential cellular interactions at the individual cellular level in ECCA, we aimed to provide a new perspective for further immunological studies and immunotherapy of ECCA. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-022-03424-5.


Introduction
Cholangiocarcinoma (CCA) is a malignant tumor originating from the epithelial cells of the biliary tree [1]. The incidence of CCA has steadily increased worldwide, and it currently ranks as the second most common hepatobiliary tumor [2], but the treatment of CCA treatment has recently shown embraced limited progression. It is widely known that CCA is characterized by a dense and reactive desmoplastic stroma around the cancer mass. Immunocytes are important components of the stroma and are closely related to prognosis [3]. Immunotherapies targeting immune checkpoint inhibition Page 2 of 15 Xu et al. Journal of Translational Medicine (2022) 20:210 have achieved breakthroughs in melanoma and nonsmall-cell lung cancer. Checkpoint inhibitors have been reported for a small cohort of CCA patients but did not show the expected success [4]. Thus, a systematic investigation of the phenotype and function of immunocytes in CCA could help to map the immune environment and provide new strategies and targets for immunotherapy of cholangiocarcinoma.
Even traditional bulk-phase sequencing has revealed that driver mutations and aberrant regulation promote progression [5][6][7], and refined identification and functional studies in the field of immunocyte subclusters in the tumor microenvironment (TME) are still lacking. The application of single-cell RNA sequencing (scRNAseq) to analyze the TME could illuminate the intricate immunocyte composition of the TME and the corresponding phenotypes and functions, particularly for lowabundance cell subpopulations. Recently, scRNA-seq has achieved breakthroughs in the study of various tumors, including lung cancer [8,9], hepatocellular carcinoma [10,11] and colorectal cancer [12]. Analysis of subclusters of mesenchymal cells at the single-cell level revealed the potential physiological functions of each subcluster and identified potential immunotherapy targets [10,13]. Currently, there are few scRNA-seq data available for cholangiocarcinoma, and we could only find a report on intrahepatic cholangiocarcinoma [14]. However, scRNAseq has never been used to study extrahepatic cholangiocarcinoma (ECCA).
Based on scRNA-seq, we captured data on 13526 single-cell transcriptomes to characterize the immune microenvironment of ECCA from the matched tumor tissues, paratumor tissues and peripheral blood of two ECCA patients. Our study clustered immunocytes, revealed immunocyte heterogeneity and described the transcriptome characteristics of different immunocyte clusters. We found potential immunotherapy targets, for example, ACP5, MAGEH1, TNFRSF9 and CCR8 for tumor infiltrating Tregs and MT1 for exhausted CD8+ T cells (CD8+ Tex cells). We also discovered some new subclusters, such as "intermediate" CD8+ Tex and "nonclassic" plasmacyte subclusters (CD27+ CD138+ CD38−). Moreover, we elucidated a potential cell-to-cell interaction network among Tregs, CD8+ Tex cells and M2 macrophages in ECCA.

Patients and clinical sample collection
Two ECCA patients who did not receive chemotherapy or radiation were enrolled in this study. Patients signed informed consent forms independently and agreed to donate tumor tissues and matched paratumor tissues and peripheral blood for scientific study. Our study was approved by the Tongji Hospital Research Ethics Committee and the Institutional Review Board. The demographic characteristics and clinical information of the patients researched by scRNA-seq are listed in Additional file 3: Table S1. The samples from two patients were collected at different times but by the same methods. After collection during surgery for ECCA, tumor and paratumor tissues were immediately transferred into MACS Tissue Storage Solution (Miltenyi Biotec GmbH), and peripheral blood was obtained before surgery. All fresh samples were stored on ice during transportation.

Isolation of mononuclear cells
Tumor tissues and matched paratumor tissues were dissociated into single cells with a Human Tumor Dissociation Kit (Miltenyi Biotec GmbH) under the guidance of the kit instructions and resuspended in HBSS (Thermo Fisher Scientific) plus 0.04% bovine serum albumin (Sigma-Aldrich). Then, blood, tumor tissue and paratumor tissue, simultaneously, were processed by the same methods. Percoll separation solutions (Sigma-Aldrich) were adjusted to isotonicity and densities of 1.077 g/ml and 1.050 g/ml and were placed successively in tubes. The single-cell suspension obtained above and the matched blood mixed with an equal volume of PBS were laid on top of Percoll and centrifuged at 600×g for 30 min at 20 °C. Mononuclear cells, mainly lymphocytes and myeloid immunocytes, were collected from the third layer between two Percoll layers of different densities. Viability was confirmed to be> 85% in all samples with an AO/PI double staining kit (Thermo Fisher Scientific).

ScRNA-seq library construction and sequencing
ScRNA-seq library construction was performed using the Chromium Single Cell 5′ Library & Gel Bead Kit (10× Genomics) for all samples according to the manufacturer's instructions. Briefly, a single-cell suspension (containing 10000 cells), barcoded gel beads and partitioning oil were loaded onto Chromium Chip A to generate a single-cell gel bead-in-emulsion (GEM). Cells were captured in GEMs during targeted cell recovery. After the reverse transcription step, barcoded cDNA was purified with Dynabeads, and PCR amplification was performed as follows. A 5′ gene expression library was constructed on the basis of amplified cDNA. For gene expression library construction, 50 ng of amplified cDNA was fragmented and end-repaired, 150-bp paired-end reads and raw BCL files were produced after the cDNA was double size-selected with SPRIselect beads and sequenced on a NovaSeq platform (Illumina). Raw data in FASTQ format were assembled from the raw BCL files using Illumina's bcl2fastq converter.

ScRNA-seq data processing and sample aggregation
We used Cell Ranger Software Suite (v.3.1.0) to manipulate the FASTQ files, aligned the sequencing reads to the GRCh38 reference transcriptome build using STAR [15] and generated a filtered UMI expression profile for each single cell. Cell barcodes were then automatically determined based on the distribution of UMI counts. The following inclusion criteria were then applied to each cell: gene number between 500 and 6000, UMI count> 1000 and mitochondrial gene percentage < 10%. Genes expressed in less than three cells were excluded. DoubletDecon was used to detect and remove doublets. After filtration, a total of 13,526 cells were retained for the following analysis. Finally, "FindIntegrationAnchors" and "IntegrateData" in Seurat v.3 were used to integrate the filtered gene-barcode matrix of all samples to remove batch effects across different donors [16]. "ccRemover" was used to identify and remove cellcycle effect [17].

Cell clustering and annotation
"LogNormalData" in Seurat v.3 was used to normalize the filtered gene-barcode matrix. Then, the top 2000 highly variable genes were calculated with the "Find Variable Features" function in Seurat. Principal component analysis (PCA) was conducted by using the top 2000 variable genes. Then, the top 50 principal components were used to perform UMAP to visualize the cells. In addition, cell clustering was carried out on the basis of PCA-reduced data for clustering analysis with Seurat v.3. Cell clusters were annotated according to SingleR [18] and marker genes; for example, the phenotypes CD3E+ CD4+ , CD3E+ CD8A+ and CD79A+ were used to identify CD4+ T, CD8+ T and B cells, respectively.

Reclustering of CD4+ T cells, CD8+ T cells, B cells and myeloid immunocytes
CD4+ T cells, CD8+ T cells, B cells and myeloid immunocytes annotated according to SingleR and marker genes were reintegrated and reclustered by Seurat v.3. Reclustered CD4+ T cells, CD8+ T cells, B cells and myeloid immunocyte subclusters were named by cluster number and cell type; for example, cluster 3 in CD8+ T cells was named C3-CD8. Specifically. With the help of canonical correlation analysis and PCA of the first 50 dimensions, CD4+ T cells, CD8+ T cells, B cells and myeloid immunocytes of all samples were integrated. The parameter k.filter was set to 120. In the clustering step, the parameter resolution was set to 0.8. The "Monocle 3" in Seurat ordered CD8+ T single cells in an unsupervised manner and was used to analyze single-cell trajectories along pseudotime to discover the developmental transitions of CD8+ T subclusters.

GO enrichment analysis
GO enrichment analysis of differentially expressed gene sets was carried out in the GOseq R and KOBAS 3.0 packages. GO terms with adjusted P values below 0.05 were considered significantly enriched among differentially expressed genes.

Cellular interaction analysis
CellPhoneDB (v2.1.7) [19] was used to classify cellular interactions for cells derived from ECCA tumor and matched paratumor tissues. Receptor-ligand crosstalk between two cell populations was identified according to the expression of receptors of one cell type and ligands of the other cell type. The expression frequency of a given receptor ligand gene was required to be > 10% in the corresponding cluster. The P value indicated the likelihood of a given receptor-ligand complex by calculating the proportion of means equal to or higher than the actual mean. The significant mean shows the total mean of the average expression values of each partner in the corresponding cell type interaction pair.
Tregs are characterized as inhibitors of antitumor immunity and promote tumor development and progression. It is of great importance to see a strong enrichment of Tregs (C5-CD4) in ECCA tumor and paratumor tissues in our study (Fig. 2C, D). Furthermore, compared with conventional T cells (Tcon, from C1-CD4 to C4-CD4), Tregs (C5-CD4) exclusively expressed 333 genes in this study, including marker genes (FOXP3   Fig. S1B).
Sixty-three genes were highly expressed in nonclassic PC subcluster C3-B, and IGHG4, IGHG1, MZB1, SDC1, DERL3, XBP1, PRDX4, FKBP11, SSR4 and IGLV3-1 were the top 10 highly expressed genes. GO enrichment analysis revealed that highly expressed genes were enriched in the pathways of immunoglobulin complex and complement activation (Fig. 4E). Analysis of C5-B showed a total of 431 highly expressed genes, and TRIB1, AQP3, B9D1, CYTOR, CNKSR1, FNDC3B, IGF1, TNFRSF17, TXNDC5 and SLAMF7 were the top 10 highly expressed genes. GO enrichment analysis revealed that highly expressed genes were enriched in the immunoglobulin complex and response to protein folding (Fig. 4F).

Discussion
It is widely known that CCA has a poor prognosis and is a challenging malignancy characterized by late diagnosis, low resectability, early recurrence, drug resistance and especially poor prognosis. Immunocytes exhibit obvious heterogeneity and promote or suppress cancer progression. Previous studies have identified immunotherapeutic targets, including PD1, PDL1 and CTLA4, that promote the treatment of non-small-cell lung cancer, melanoma and renal carcinoma. Deeper insight into the heterogeneity of immunocytes in CCA is urging, and it is promising to provide potential immunotherapeutic targets. Herein, based on scRNA-seq, we analyzed immunocytes of ECCA patients to describe the transcriptome characteristics and potential function of each immunocyte subcluster. As CD4+ Treg cells promote tumor immunosuppression and are associated with poor prognosis, novel immunotherapy targets for Treg cells may improve the treatment of ECCA. Our study analyzed Treg cells and revealed 333 exclusively expressed genes. Eighteen genes identified in our study were consistent with those reported by earlier studies in melanoma, breast, colon, lung and liver cancers, including coinhibitory receptors CTLA4 and TIGIT and costimulatory receptors TNFRSF4 and TNFRSF18. Previous studies have proven that CTLA4, TIGIT, TNFRSF4 and TNFRSF18 cannot promote and instead suppress the proliferation and immunosuppression of Tregs [45][46][47]. Among the 18 Treg-specific overlapping genes, ACP5, MAGEH1, TNFRSF9 and CCR8 were exclusively expressed in tumor-infiltrating Tregs and are potential immunotherapy targets for tumor-infiltrating Tregs.
In our study, CD8+ Teff subclusters in peripheral blood, not in tumor/paratumor tissues, highly express CX3CR1 which is related to strong cytotoxic effector function, unique migration pattern and positioning adjacent to the entry site of pathogens [29]. Interestingly, the expression levels of the cytotoxic effectors PRF1, GNLY, GZMA, GZMB and GZMH in CD8+ Tex subclusters were even richer than those in CD8+ Teff subclusters, which implied that CD8+ Tex cells still tended to respond to cancer cells. However, the CD8+ Tex subclusters, characterized by high expression of coinhibitory receptors, presented limited proliferation and differentiation, which suggested an immunosuppressive TME in ECCA. Moreover, our study revealed two CD8+ Tex subclusters (C1-CD8 and C4-CD8). Compared with subcluster C4-CD8, C1-CD8 weakly expressed coinhibitory receptors but highly expressed MT1 which correlates with loss of effector function and acquisition of a dysfunctional phenotype of CD8+ T cells [28]. In addition, pseudotime trajectory analysis confirmed that C1-CD8 could differentiate into C4-CD8. Thus, C1-CD8 may consist of "intermediate" CD8+ Tex cells and could differentiate into "complete" CD8+ Tex cells C4-CD8.
Consensus exists the utility of combined assessment of CD38 and CD138 for the identification of plasma cells, and CD38 emerged as a highly valuable PC-identification marker and can be reliably used for the identification of both normal and myeloma PCs. In our study, in addition to the classical PC subcluster coexpressed marker genes CD138 and CD38, a rarely described "nonclassic" PC subcluster with low levels of CD38 was observed. GO enrichment analysis revealed that the "nonclassic" PC subcluster was associated with the pathways of immunoglobulin complex and complement activation, which may be related to the regulation of innate immunity.
M2 macrophages are involved in the anti-inflammatory immune response and promote tumor activity. Compared with peripheral blood M2 macrophages, tissues/paratumor tissue M2 macrophages highly expressed the marker gene MCR1 (CD206) and the protumorigenic genes TGFβ, IL10, VEGFA and MMP9, which indicated stronger protumorigenic and immunosuppressive properties, and MCR1 may be a better marker to filter tumor-infiltrating M2 macrophages. pDCs are myeloid immunocytes but show both lymphoid and myeloid lineage transcriptome characteristics. In addition, pDCs highly express SCT and may have endocrine functions. pDCs may form an unclear cluster and need to be further studied.
A potential cell-to-cell interaction network among Tregs (C5-CD4), CD8+ Tex cells (C1-CD8 and C4-CD8) and M2 macrophages (C6-Mye) in ECCA was revealed in this study, implying the formation of a suppressive immune microenvironment in tumors via potential communication among diverse immune clusters. Consistent with previous reports, CCR8 is hyperexpressed on Tregs, and its ligand, CCL4, secreted by M2 macrophages, has been demonstrated to recruit Tregs to the tumor niche in other cancers [48,49]. Moreover, connections between TIGIT-NECTIN2 and CTLA4-CD86 have been reported in previous studies, by which Tregs could regulate their own activation [50]. HLA-E-KLRK1 was predicted to engage in cross-talk between Tregs and CD8+ Tex cells, by which CD8+ Tex cells could escape being neutral and promote tumor progression [38]. Notably, CCL5-CCR5 and CCL4-CCR5 were predicted to make great contributions to interactions among Tregs(C5-CD4), CD8+ Tex cells (C1-CD8 and C4-CD8) and M2 macrophages(C6-Mye), which were characterized by recruiting and promoting activity of tumor associated immune cells, including TAM, Th17 and Treg cells [51][52][53]. In addition, receptor-ligand binding involving VEGF is densely enriched in monocyte clusters, which could contribute to tumor angiogenesis. All the potential links between diverse types of immunocytes are crucial for the preservation of homeostasis and functional regulation in the TME. Immunosuppressants, such as PD-L1 and CTLA4 inhibitors, have been moved into clinical trials for ECCA treatments, but they have not met expectations [54]. Based on the receptor-ligand analysis in this study, the immunosuppressive cross-talk of HLA-E-KLRK1 and LGALS9-CD44 between Tregs and CD8+ Tex cells suggests that they may become potential immunotherapeutic targets in ECCA. There are several limitations to the present study. First, the sample size in our study is quite small. Compared with some solid tumors well and much studied by scRNA-seq, including intrahepatic cholangiocarcinoma, hepatocellular carcinoma, breast cancer and lung cancer, the samples, particularly paratumor tissues (weigh < 0.1 g), which we can acquire from extrahepatic cholangiocarcinoma (ECCA) patients are too small to isolate enough qualified mononuclear immune cells, as the high sample quality (single cell count > 10000 and viability > 85%) was needed to perform 10X Genomics scRNA-seq. Besides, functional assays are needed to prove our results at the protein level. However, we first clustered immunocytes and revealed immunocyte heterogeneity in ECCA-matched tumor tissues, paratumor tissues and peripheral blood based on scRNA-seq. We found potential immunotherapy targets, for example, ACP5 and CCR8 for tumor-infiltrating Tregs and MT1 for CD8+ Tex cells. We also described some interesting subclusters: "intermediate" CD8+ Tex cells and the "nonclassic" PC subcluster. In addition, we proved the existence of a potential cell-to-cell interaction network among Tregs, CD8+ Tex cells and M2 macrophages in ECCA. Our study may provide a new perspective and potential targets for the treatment of ECCA.